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ABSTRACT 

We present an efficient method for two dimensional inversions for tfie solar rota- 
tion rate using the Subtractive Optimally Locahzed Averages (SOLA) method and 
a modification of the ® technique proposed by Sekii (1993a,b). The SOLA 
method is based on explicit construction of averaging kernels similar to the Backus- 
Gilbert method. The versatility and reliability of the SOLA method in reproducing 
a target form for the averaging kernel, in combination with the idea of the ® M 
decomposition, results in a computationally very efficient inversion algorithm. This is 
particularly important for full 2-D inversions of helioseismic data in which the number 
of modes runs into at least tens of thousands. 

Key words: Sun : oscillations of - Sun : rotation of - Sun : structure of - Numerical 
methods 



1 INTRODUCTION 

The solar 5-minute oscillations can be described as a super- 
position of eigenmodes of non-radial pulsation. Each mode 
is identified by three integers (n, Z,m), where I and m are 
the degree and order respectively of a spherical harmonic, 
and n is essentially the number of radial nodes in the dis- 
placement eigenfunction: m can take all values from —I to 
In a spherically symmetric, nonrotating star, the fre- 
quency of an eigenmode would be independent of m and 
thus there would be multiplets of {21 -\- 1) modes with iden- 
tical frequencies, each multiplet corresponding to an (n, I) 
pair. Rotation lifts this [21 + l)-fold degeneracy. The dif- 
ference in frequency between modes in the same multiplet 
is called the (rotational) frequency splitting. The frequency 
splitting is determined by the rotation rate inside the Sun 
and can be used in an inverse problem to probe the Sun's 
internal rotation. In particular, they enable one to perform 
2-D inversions for the rotation rate as a function of radius 
and latitude. 

Large helioseismic data sets should soon be available 
from various observational campaigns, notably the Global 
Oscillations Network Group (GONG) and MDI-SOI on 
board the SOHO satellite. To make optimal use of these 
data, the algorithms for full 2-D helioseismic inversions need 
to become efficient in handling several hundreds of thou- 
sands, or even millions, of data (modes) simultaneously. 
Optimally localized averages (OLA) techniques, which have 
proved very popular for 1-D helioseimic inversions involving 
only a few thousand data, require a matrix to be inverted 



whose order is the total number of data. This is prohibitively 
expensive computationally in the 2-D case - though one may 
be able to make the computation tractable by preprocessing 
to reduce the number of data to which the OLA is applied 
(Christensen-Dalsgaard et al. 1994). Least-squares tech- 
niques (e.g. Schou et al. 1994), which require a matrix to 
be inverted whose order is the number of base functions, 
are also expensive in the 2-D case, where one might want a 
discretization of e.g. 200 bins in radius and 100 in latitude 
which gives a matrix of order 2 x 10*. 

Sekii (1993a,b) has exploited the fact that the rotation 
kernels are nearly separable in radius r and colatitude 9 to 
develop a so- called Wt} » Wt} inversion technique in which 
the true kernels are approximated by ones which are exactly 
separable. This results in a problem where the order of the 
largest matrix one has to invert is only the number of (n, l)- 
multiplets, which is only a few thousand. 

Here we propose a modification to Sekii's ® IR^ 
method, in which the small deviations from separability of 
the kernels are taken into account. The computational bur- 
den is the same as for Sekii's approach; hence it is just as 
efficient. Again the problem is reduced to a series of ID 
inversions, for which we use the subtractive optimally local- 
ized averages (SOLA) method of Pijpers & Thompson (1992, 
1994; hereafter PTl, PT2). The SOLA has the advantage 
over other ID inversion methods that it is possibly to keep 
close control over the averaging kernels that it produces: 
this turns out to be important for the ® inversion, 
as we shall see below (cf. also Sekii 1994). 
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2 THE ig) JR^ METHOD 

In a spherically symmetric, nonrotating star, the frequency 
LOni of a spheroidal mode of oscillation of radial order n and 
degree I is independent of the mode's azimuthal order m, 
and the displacement eigenfunction is 



d L ^rini d 



d<j> 



Pr (cos 9)e' 



(1) 



with respect to spherical polar coordinates {r,6,(f)): here 
^ni and rini are calculable functions, given a solar model; 
L = + 1); and is the Legendre function of degree I 
and order m. 

For a slowly rotating star the frequencies also depend 
on m. According to linear perturbation theory, which is an 
excellent approximation for the Sun, the difference between 
frequencies of modes with the same values of n and I, but 
opposite m, is given in terms of the eigenfunctions of the 
nonrotating star by 



where 



Dr, 



2mD„im 



Kr,im{r,e)n{r,e)drdcose. 



(2) 



Here and in the following, we have made the radial variable 
r dimensionless by dividing it by the surface radius. 

An OLA inversion for the rotation amounts to finding 

coefficients {c„i,n(ro, 9o)} so that 



n(ro,6'o) = Cnlmjro, 9o)Dnlm = 



j j CnlmKnlm^ dr d COS ^ ■ 



(3) 



is a localized average of the actual rotation rate f2 near r = 
ro,0 = 00. To find the coefficients with a naive application 
of OLA would require a matrix to be inverted whose size was 
the total number of observed eigenmodes i.e. all available 
{ri,l,m) combinations. This is prohibitively expensive with 
a very laxge mode set. This paper is concerned with finding 
suitable coefficients in a computationally less expensive way, 
by exploiting properties of the kernels Knim- 
The Knim arc given by 

Kr.imir,e) = Fr\r)G'rie) + F^\r)Gir{e) , (4) 

where 

Fi"' = piry[eniir)-2L-'Uir)Vniir) + vUr)]/Ini (5) 

= p{ry[vUr)]/Ini, (6) 

p being the density and and £,„i , ri„i the components of the 
displacement eigenfunction in the nonrotating star, and 



-I 



<ir pr^i(,li+vli) ■■. 



and 



(7) 



(8) 



L-(l. 



(9) 



Note that our definition of rjni differs from that used by Sekii 
(1993a) by a factor of L. The ratio of the amplitudes of our 
^ni and rini is roughly the ratio of the vertical wavenumber 
to the horizontal wavenumber (since for p modes the waves 
are almost longitudinal): thus for p modes, which propagate 
more nearly vertically than horizontally except near their 
turning points, is larger in magnitude than rjni. More 
quantitatively, using the simplest asymptotics, the ratio of 
the amplitudes of F"' to F2' is of order uj^r^/L^cs in the 
acoustic cavity of the mode, oj being the mode's frequency 
and cs the local adiabatic sound speed. Although this is of 
order unity near the lower turning point, it is much larger 
than unity elsewhere in the acoustic cavity. The angular 
functions G'^, G^^ 

are of similar magnitude, since the two 

derivatives of Gi"" each produce a factor of L and these 
cancel the factor of L^^ . Hence, as Sekii (1993a, b) observed, 
the contribution from F^^G'^J^ is generally small compared 
to that from Fi"'G'i"', and so the original problem can be 
approximated by the separable problem 



2m 



Fi\r)G'i" (u) drdu . 



(10) 



A further aspect is that G{^ is everywhere positive, whereas 
G2" oscillates to both positive and negative values - increas- 
ingly so for larger values of {I — \m\) and this too helps 
make the integrated contribution from the F2^G2" small in 
general. 

The essence of Sekii's idea for exploiting the separability 
is as follows. (We shall work within the framework of SOLA, 

which is the ID method we shall employ.) Firstly, for each 
I one seeks coefficients (do) such that 



j2crwG'r{e) « Ti{u -uo).^ 



(11) 



where wo = cos^o, and Ti{u — uo) is some chosen target 
form that is peaked around 6 = 9o and small elsewhere. 
(The means by which such coefficients may be sought is 
described in PTl, PT2.) Then 

I 

= J I Fr' |^cr(eo)C?'r(0) j fl{r,e)drdu (^2) 



Fr'{r){Q}l''°>{r)dr, 



where (n)*""' ^ / Ti{u - uo)n{r,e)du. The second and 
final step is then to choose further coefficients c"'(ro) such 
that 



^c"'(ro)Fr'(r) « T{r-ro) 



(13) 



where similarly T(r — ro) is a chosen target function that is 
localized about some radius, ro. Now if Ti(u — wo) were in 
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fact independent of /, so (O); = (Q) say, eqs. (12) and (13) 
would imply that 



^c"'(ro)cr(^o)i?ni^« / ^c"\ro)Ff{r 



{n)'''«Hr)dr 



(14) 



where {{0))^''°'"°^ ~ // ^('^' - ro)T{u - uo)Vl{r,e) ArAO. 
This then completes the R,^ ® inversion. 

There are some drawbacks to the simple ® pro- 
cedure outlined above. One is that the linear combinations 
on the loft-hand side of (11) need to be essentially indepen- 
dent of /. A practical matter is that choosing the angular 
target functions to be independent of I does not guarantee 
that X]cr(6'o)G'r(6') will itself be independent of the de- 
gree. However, SOLA is better than other commonly-used 
methods in forcing the linear combination to accurately re- 
semble a given form (Dziembowski et al. 1994, Sekii 1994). 
A second point is that, in making all the target functions 
independent of I and ro, the angular resolution implied by 
(11) is the same at all target radii. But because we have 
many more m values for the shallowly penetrating high- 
degree modes than for the deeply penetrating modes of low 
degree, we should be able to achieve much better angular 
resolution in the outer layers of the Sun than in its deep in- 
terior. [This point was appreciated by Sekii (1993b); but one 
of the contributions of the present paper is to suggest how 
the angular resolution might be chosen appropriate to the 
target depth.] Thus we need to relax the restriction that the 
combinations of latitudinal kernels are independent of both 
degree and radial target. In order to preserve the property 
that the radial inversion can be treated as a 1-D inversion, 
we shall allow the target functions Ti, and hence the coeffi- 
cients q"*, to depend on the radial as well as the latitudinal 
location of the target point. Thirdly, although it may well 
be an excellent approximation to neglect the F;^'G2" terms 
for most presently observed modes, we hope that forthcom- 
ing datasets will extend to even lower frequency low-degree 
modes, for which Sekii's approximation is less good. For this 
reason, we wish to include the F-^^G^J" contribution to the 
kernels while retaining the advantage of the basic R^ (8 R'^ 
concept. 

We therefore now introduce a modification to Sekii's 
method to take into account the terms when trying 

to form a localized average. Specifically, in the latitudinal 
localization, we seek coefficients Ci^{ro,9o) to minimize 



+ 2^ E^ wCi C( . 



(15) 



Here 



Mo 



/ J m m' 



(16) 



Ho being an adjustable parameter; E^'^, is a suitable error- 
covariance matrix (see below). Minimizing the error term 
and obtaining a well-localized kernel are opposing aims (e.g. 

CDST) and the error weighting parameters fio arc used to 
obtain a compromise between them. We simultaneously 
force the linear combinations of Gi" and of G2^ to resemble 
the chosen target functions Tqi and , which depend on 
I, ro and ^0; and f3i is an adjustable parameter that weights 
the relative importance of matching Tg^ and matching Tq^- 
It is convenient to define function vectors 



and 



Tg2 



(17) 



(18) 



Note that we shall generally omit the superscript I on Q and 
Tg, except when we wish to emphasize their dependence 

on the degree I. 
suggestive form 



Then expression (15) can be written in 



1 

/ 



du{g - TGf[ I ^° ]{G - Tg) 



(19) 



(/) ~'m~m' 
m m' 



In effect, we are seeking to localize a vector of functions 
rather than a single function, so that the inversion can once 
again be carried out as a sequence of two 1-D inversions but 
without the approximation Sekii (1993a,b) employed. 

If a clean two-dimensional localization is finally to be 
achieved, it is desirable that the components Tc^ and Tg2 
of Tg should have similar shapes. Thus we choose to make 
the second component the same as the first, up to a multi- 
plicative scalar. 

For convenience later we also impose an exact constraint 
on the first component of G, that it be normalized such that 



1 

/ 



dug[{eo,e) 



If the latitudinal localization is successful, then 



(20) 



Ycriro,eo)Dnim = J J F"' ■gn{r,e)drdu 



^11 



F"' • Tc^ir, 6) drdu 



{Fr'+C'F^')TG,ndrdu, 



(21) 



where C"' is essentially the ratio of to Tgi • 

The second and final step in the two dimensional inver- 
sion is to find a second set of coefBcients such that : 



nl 



c"'(ro) [Fr'{r) + C'F^'{r)] 



(22) 



is peaked around r = ro and is small everywhere else. 

Actually the ("^ that enters in (22), which is defined 
explicitly below in equation (29), is not precisely the ratio 
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of the components of the latitudinal target function we use, 
though they are closely related. For our latitudinal localiza- 
tion we use 



\ feAtL 



/ u - up \ 



\ 



(23) 



make the total integral of Tqi equal to unity. 
The radial target function is 



Tf = 



/r-A, 



exp 



The expression to minimize is now 



1 



nl n'l' 



(24) 



(25) 



The free parameters A^, Ag should in general be func- 
tions of target position. The choice of A,., Ae at different 
radii and at different latitudes is an important issue. Wc 
have used the natural scaling that was presented in PT2 (cf. 
Thompson 1993) for radial and latitudinal resolution: 



A^(ro) 



8 



cs(ro) 



Ar 



(26) 



Ae(ro,6'o) =a0 — ^/l + t-ul 
ro 

where Ur and ae are constants of proportionality, indepen- 
dent of radius. To avoid problems at mq = 1 we include the 
small number e = 0.075. The strategy is to optimize the tar- 
get widths at one target location and then use the relation 
(26) to calculate the target widths for all other locations. 
One should note that it is at this point that the dependence 
of the angular resolution on radius is introduced. The choice 
of functional dependence in (26) expresses the fact that one 
expects the attainable physical horizontal resolution and ra- 
dial resolution to be very similar at all radii. 

The factor that appears in the radial inversions (22) 
is necessary because the values of tfie two latitudinal aver- 
aging kernels at the location of their maxima differs due to 
the relation (9). The factor can compensate for this in 
the radial inversions so that in the final reconstructions the 

do not suddenly acquire a much smaller or larger weight 
due to the multiplication by the G2™. Using (9) and (23) 
the ratio of the values of the two latitudinal kernels at their 
maximum is approximately : 



32(60,9 = 00) 

gi{9o,e = eo) 



(i-^g) 

A2L2 



(27) 



The dependence of C"' on the latitudinal resolution width 

implies that a straightforward application of this recipe 
means that the radial kernels to be used in the inversion 
depend on the choice of the target localization radius. This 
would defeat the purpose of using SOLA, because it would 
require calculating and inverting a new matrix for each new 
target position ro. However, one can make use of the fact 
that any mode will be used primarily at target radii ro close 
to the turning point of that mode r*. In other words the 
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Figure 1 Flowchart showing the principal steps in the modified 

(B> method presented here. 

coefficient for this mode is only large if the averaging ker- 
nel to be constructed peaks near the turning point of that 

mode. Using this information and the scaling law for the 



latitudinal resolution widths 
92(60,9 = 9o) 



(^"' becomes 



gi(6o,9 = 9o) 



/^^iwy ^28) 
\naear v J 



If the absolute value of this factor is allowed to become much 
larger than unity it turns out to over-compensate for the 
effect of the amplitude of Q2 in the radial inversion. Essen- 
tially what we wish to do with (^"' is to malce sure that the 
product F2G2 will indeed be treated as a small correction 
term, and hence should be less than unity. Thus in prac- 
tice we achieve this while retaining the (^-dependence in (28) 
by defining 



\ V ) 



(29) 



Equation (29) is used for all except the lowest degree modes 
for which the turning point is not close to the point at which 
the mode is actually used in localized kernels. 

The factors Pi are used to compensate for the difference 
in the absolute value of the integrals of Tg-^ and Tg^ . We 
take Pi = (Agl/^) where the () denotes a simple average over 
the latitudinal target points. 
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Now combining (21) with (22), this time keeping track 

of data errors e„im, yields : 



c"'(ro) ^ cr(ro, eo)Dn 



1 

y dr ^ c"' [F^\r) + C'F^\r)] {n(r,eo))' + 



rtl m 



1 

= / dr^(ro,r){n(r,eo)) + ^c"'^cre„,. 



rtl m 



The constraint 
1 

J drJ^{ro,r) = 1 



(30) 



(31) 



is now sufficient to ensure that the complete 2-D averaging 
kernel that is constructed is correctly normalized : 

1 1 

dr J duK:{ro,eo;r,e) = 1 . (32) 
-1 
where 

nl m 

(33) 

The difference between the target form and the actually con- 
structed kernel will give rise to a 'systematic' error a measure 
of which is x ■ 



1 1 
X = J drj du 
-1 



\ \ nl — m i jpnl/^lm , T-ini v^imN 



-TpTai 



The error covariance matrix E^^^, used in (15) and (19) 
is derived from that of the observed frequencies. Since the 
localization in colatitude is done at fixed I, the indices m and 
m' run over the range of m, viz from 1 to L The appropriate 
covariance matrix is therefore Erii,n,.nim' where / is the degree 
under consideration and n can be any order such that the 
(n, I) multiplet is in the dataset. The overall magnitude of 
E is irrelevant, since this is taken out by the scaling in /i. 
One should note that the error covariance of m and m' for 
a given I is taken to be independent of the value of n. This 
restricts the permitted data error covariance matrices. Of 
course one can proceed even if the errors do not satisfy this 



(34) 




0.2 0.4 0.6 



Figure 2 Constructed unimodular averaging kernels for various 
values of latitudinal and radial resolution widths. Each panel is 
labelled with the associated error in the solution, in nHz. Middle 
panel ar = 8.0, ag = 1.5. Upper left ar = 6.0, = 1.5, up- 
per right Qr = 8.0, ae = 1.0, lower left ar = 10.0, = 1.5, 
lower right ar = 8.0, ag = 2.0. The target location is 
ro = 0.7-R, ^0 = 45° in all cases. The highest contour is 
~ 10% below the peak value and the other contours are spaced 
in 1/8 of that value. The crossbar has an extension that corre- 
sponds that of a Gaussian at the lowest contour level which is 
at |r - ro| = 1.478Ar, \u - uo\ = 1.478A(j. All error weighting 
parameters no =0.1 

assumption, but the inversion will be less "optimal" . 

The error covariance matrix E„( „/;/ in the radial inver- 
sion is also derived from that of the observed frequencies. 
We assume that the data error covariance matrix can be fac- 
torised, into an error (co-)variance between multiplets and 
an error (co-)variance for m given any fixed I. If this is not 
in fact the case, one can still approximate the errors in this 
manner, but once again the inversion will presumably not 
be optimal. In the radial inversions the matrix is 
obtained by : 

ve(" ,crcr X 

/ V m m' ^ I 

(35) 



nl n'l' 



1/2 



Here eni„'i' is a matrix which expresses the error variance/ 
co-variance for all (n, i)-multiplets. We conclude this pre- 
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Figure 3 Unimodular averaging kernels at radii ro = 0.1,0.4,0.7,0.95 and co-latitudes 0o = 15°, 30°, 60°, 90°. Target locations are 
indicated by crosses; standard errors are as quoted. Localization parameters were ar = 8.0 ag = 1.5 and error weighting factor no = 0.1 
for all kernels. 



sentation of the modified (g M} method by outlining the 
principle stops in a flow chart (Fig. 1). The flrst and second 
box represent stops that aro preparatory to the actual in- 
version, the third box represents the latitudinal part of the 
inversion and the fourth box the radial part. The final box 
represents the process of combining the linear coefficients 
from the inversion with the data. 



3 RESULTS 

To illustrate the results which may bo obtained with the 
method, we have applied it to artificial splitting data. The 
modeset is exactly the same as used in the GONG Hare and 



Hounds exercise (Gough & Toomre 1993): in brief, the set 
consists of f and p modes below 5mHz with 

/ = 1, 2, . . . , 16, 18, . . . , 50, 55, . . . , 150, 160, . . . , 250 

making a total of 69662 individual modes with positive m 
values (in 1380 n,l multiplets). The assumed uncertainties 
on the frequency differences ujnim — (^ani were of the form 
ai{i') = f{y)g{l), the functions / and g being estimated from 
Fig. 3.5 of Gough & Toomre (1993), and so are roughly the 
same as the assumed uncertainties in the Hare and Hounds 
exercise. The range of ai{v) was from about 4nHz at low 
degree and frequency to roughly 700nHz at high degree and 
frequency. Note that under these assumptions, the uncer- 
tainties in the i?nim are f{i>ni)g{l)/m. 
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Figure 4 Unimodular averaging kernels at radii ro = 0.1,0.4,0.7,0.95 and co-latitudes 9o = 15°, 30°, 60°, 90°, without also localizing 
the F2G2 -terms. Target locations are indicated by crosses; standard errors are as quoted. Trade-off parameter values are the same as 
for Fig. 3 



It should bo noted that the GONG Hare and Hounds 
dataset contains some low frequency p and f modes for which 
the assumptions about the location of the lower taring point 
made in the previous section are inaccurate or inappropriate. 
Specifically, for f modes with I < 13 and and also for the 
I — l,n = 1 p mode, is computed using (27) and (26) 
with To set to 0.1 instead of the turning point radius rt- 

The resolution and localization achieved in the in- 
versions can be seen by inspecting the averaging kernels 
lC{ro,6o;r,9), defined in equation (33). In Fig. 2 we illus- 
trate averaging kernels at ro — 0.7R, 9o = 45°, for different 
choices of target widths Ar and A^. As can be judged from 
the crosses superimposed on the kernels, the constructed 



averaging kernels match closely the specified target forms. 
The centre kernel has a radial width Ar of 0.0648 (frac- 
tional) solar radii and an angular width Ag of 7.75°. This 
yields a standard error in the inversion of 3.3 nHz. Compar- 
ing this level of uncertainty with the range of the rotation 
rate observed at the Sun's photosphere (roughly 320 nHz - 
460 nHz), or the range (approximately 300 nHz - 1500 nHz) 
of the rotation rate present in our artificial example (see Fig. 
6a), it is evident that this is an acceptably small error for 
many purposes. The resolution could be squeezed further 
in cither radius or latitude, but at the expense of increasing 
the standard error. Conversely, the error could be reduced 
somewhat by degrading the radial or latitudinal resolution. 



8 F.P. Pijpers and M.J. Thompson 



A more complete picture of what may be achieved glob- 
ally is provided by Fig. 3, which shows averaging kernels at 
various target radii and latitudes. These were all obtained 
with the same values of the scaling parameters Ur — 8.0, 
ag = 1.5, as well as the same value of the error trade-off pa- 
rameter fiQ = 0.1. As can be seen, the resolution is best in 
the outer part of the Sun, as expected from the asymptotic 
scalings (26) , because modes with high I are sensitive to the 
rotation in these layers. The standard errors arc also small- 
est in this region. The averaging kernels are very cleanly 
localized near the intended target location, for most loca- 
tions in the outer half of the Sun. Only near the pole does 
the kernel exhibit noticeable nonlocal structure - only rela- 
tively low-m modes sample the region close to the pole, so it 
is harder to construct a kernel there. The standard error is 
also much higher for the near-polar inversion, for the same 
reason. 

Even for target locations as deep as 0.4i?, the equatorial 
kernels arc very well localized, and the error is only 4nHz. 
The resolution could of course be improved further if the 
error was allowed to increase but, as illustrated in Fig. 2, 
one could not expect to improve the resolution substantially 
without at least doubhng the standard error. The polar 
kernel at this radius shows considerable nonlocal structure, 
though it still has a main peaJj close to the intended lo- 
cation. Thus it is possible to localize kernels at different 
latitudes and thus to achieve latitudinal resolution at this 
depth, albeit with a considerable uncertainty a. This is not 
the case at ro = O.li?, our deepest target radius. Here all 
the kernels are similar and are localized at low latitudes. 
Nonetheless, apart from an obvious positive region in the 
convection zone, and a thin negative region near the surface 
at high latitudes, these kernels are reaonably well-localized. 
However, as discussed below, they depend heavily on the 
low-degree f modes in the Hare and Hounds dataset, which 
have so far not been observed. 

To illustrate the importance of the F.f G'j", we have 
recomputed the coefficients used to construct the kernels in 
Fig. 3, setting f3i = and C"' = for all n and I. The values 
of all other parameters are unchanged. The results should 
still be superior to a basic IR} (g)R^ inversion, as described at 
the beginning of Section 2, because we are choosing the coef- 
ficients cJ™ as a function of ro as well as of 9o. As expected, 
we see in fig. 4 that neglecting the F^^G^" contribution to 
the mode kernels has little effect in the outer part of the Sun, 
but degrades the kernels at ro = 0.4i? at all but the equa- 
torial target. At the deepest target locations, the kernels 
constructed without taking F2^G'^ into account are very 
poor indeed. This is precisely as one would expect, because 
for such a deep target location the inversion needs to use low- 
degree (and, if they are available, low-frequency) modes, for 
which it is a poor approximation to neglect 

This point is further illustrated Fig. 5, where the radial 
coefficients c"' are shown for two target locations, using the 
same parameters as for Fig. 3. Panel (a) shows coefficients 
for a target radius ro = 0.7R: in that case, the coefficients 
are approximately a function of lower turning point, as has 
been seen in 1-D inversions (Christensen-Dalsgaard et al. 
1990). At ro = O.l-R the low-degree f modes (n = 0) and 
gravest p modes are given most weight, as can be seen in 
panel (b) . For such modes it is a poor approximation to ne- 
glect the ^^'02"" contribution to the mode kernels, which is 
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Figure 5 a. The coefficients c"' for radius ro = 0-7, 60 = 45° 
as a function of turning point, b. The coefficients c"' for 
ro = 0.1, Oo = 90° as a function of degree n 

why the kernels in Figs. 3 and 4 differ so markedly for such 
deep target locations. It could be argued that the f modes 
should be excluded from the artificial datarsets because they 
have not yet been observed. It is clear that localizing a ker- 
nel as deep as ro = 0.1 J? would then be more difficult since it 
is clear from fig. 5b that it is primarily those that arc used in 
the localization this deep. At ro = OAR the coefficients for 
the f modes are already of the same magnitude and smaller 
as all the other coefficients. Their omission should therefore 
not significantly affect the quality of the kernel nor the error 
magnification. It is clear from a comparison of figures 3 and 
4 at ro = OAR that the improvement with the new method 
proposed here is still appreciable, because of the improved 
treatment of the gravest p-modes. 

Although the averaging kernels contain complete infor- 
mation about the resolution of an inversion, it is nonetheless 
instructive to see how well our method reconstructs a par- 
ticular rotation profile. For this purpose, we have invented 
a rotation profile, shown in Fig. 6a. Its features include 
two jets, a fast-rotating core and polar region, and an en- 
hanced rotation at r ~ 0.9ii. Using the appropriate mode 
kernels, the rotation profile was then used to compute ar- 
tificial splitting data Dnim, which we then inverted using 
the same parameters as for Fig. 3. In the inversion, wc as- 
sumed the error uncertainties described at the beginning of 
this section, but we performed two test cases, one where no 
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Figure 6 a. The original artificial rotation rate. b. The reconstruction using error-free data. c. The reconstruction with noisy data. d. 
The errors in the reconstructed rotation rate of c. The parameter vaues used in the reconstructions were ar = 8.0 ag = 1.5 and no = 0.1. 



errors were actually added to the data and a second where 

independent Gaussian errors with the assumed standard de- 
viations were added as noise to the data. The inversion of 
the noise-free data is shown in Fig. 6(b). The fast-rotating 
core, the abrupt change in rotation at r w 0.7R and the 
general trend of the rotation in the outer part of the Sun 
are all correctly inferred. The outer jet is rather poorly re- 
solved, and the deeper jet is not convincingly detected at 
all. This is consistent with the width of the main peak of 
the averaging kernels shown in Fig. 3. In such a noise- 
free case, one could of course make much narrower kernels 
and hence obtain better resolution, because data errors arc 
not a concern. However, the main purpose of showing the 
noise- free inversion is so that it can be compared with the in- 



version of noisy data in Fig. 6(c). Except near the pole, the 

corresponding contours in the two panels can very clearly 
be identified. The relatively small differences are due to the 
propagation of data errors, and are consistent with the stan- 
dard errors quoted in Fig. 3 and with the contour plot of 
standard errors in Fig. 6(d). It is clear from panel (d) that 
the standard crrrors increase substantially in the near-polar 
region, which accounts for the greater distortion of the noisy 
inversion near the pole. In the rest of the Sun, however, 
one can see from comparing panels (a)-(c) that it might be 
worthwhile to choose smaller widths for the averaging ker- 
nels, thus improving the resolution, even at the expense of 
somewhat greater errors. Fig. 7 illustrates that one could 
make a statistically significant detection of the deeper jet in 
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Figure 7 The effect of improved resolution on the error and the 
deduced rotation rate at the location of the 'jet' 



this way. 

Fig. 8 illustrates the difference between the inferred 
rotation from the noisy data and the true rotation profile 
divided by the standard error. In rather few places does the 
deviation between inferred and true rotation rise above 2a. 
One place is the core, even though the inversion reproduces 
qualitatively the fast rotation. The other most significant 
deviations are around the two jets, particularly the deeper 
one. In the latter case, the broad averaging kernels have 
spread the jet out, so that not only is there a strong negative 
deviation where the prograde jet should have been, there is 
also a positive deviation around the jet where the inferred 
rotation is higher than the true rotation. Also shown in this 
figure are the locations where we constructed the inversion 
solution upon which the contour plots in Fig. 6 were based. 




r 



Figure 8 The difference between reconstructed and original, di- 
i^ided at each point by the local a. Dots mark the central locations 
)f all averaging kernels used in the reconstruction 



i DISCUSSION 

We have previously demonstrated (PTl, PT2) that the 
sOLA method produces kernels that are of as good a quality 
is those obtained with the more traditional OLA formula- 
tions and shown that the SOLA method is computationally 
more efficient. 

Our (gi TR^ nversion differs from the one proposed 
Dy Sekii (1993a,b). He approximated each mode kernel as 
% product Fr'(r) and 0^(9), where Ff', Gf are given by 
squations (5) and (8). The major difference of principle is 
that Sekii has thus approximated the kernels, whereas we 
ieep the small term _F^'G2' and therefore make no approx- 
mation (beyond assuming that the rotational splitting is 
correctly described by first-order theory). Sekii's original 
method would not be recovered completely from ours in the 
limit Pi 0, C"' ^ 0, because we do adjust latitudinal res- 
olution as a function of depth which Sekii (1993b) already 
suggested but did not yet implement. 

For all but the low degree modes (say / ^ 5), Sekii's 
approximation of neglecting F^' should be very good. How- 
ever, to perform inversions in the deep interior it is highly 
desirable to include the lowest degree modes, which pene- 
trate deeply into the star, and this is apparently a draw- 
back of Sekii's approximation (Sekii 1994). Now the effect 
on the inversion in the deep interior of neglecting the -F"' 
term is uncertain a priori, since even for low degree modes 
the component is not large compared to F"'. However, 
from figures 3 and 4 it can be seen that there is a substan- 
tial difference between kernels at ro = 0.1 depending on 
whether or not the F2G2 terms are included. At this radius, 
the present method makes substantial use of the low-degree 
n = modes (Fig 5). It could be argued that data for such 
modes are not available at present. But if they are observed 
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by e.g. GONG, they will be enormously helpful for probing 
the core region and it is highly desirable that our inversion 
methods make accurate use of them. 

Apart from retaining the F2G2 term, and taking this 
into account in the latitudinal and radial localizations using 
SOLA, the present approach is very much in the same spirit 
as Sekii's orig inal (g) method (Sekii 1993a,b). Subse- 
quently, Sekii has made a modification to the technique, to 
make a 2-D integral in the radial localization (Sekii 1993b, 
1994). Such a modification might be included in our ap- 
proach also. 

Any multiple of Fj^'Gi"' could be added into the first 
term on the right-hand side of equation (4) and maintain the 
form of the right-hand side, provided the same multiple was 
subtracted from the second term. In particular wo could 
add to the first term and subtract it from the 

second. This is equivalent to replacing with Fi"'-FC"'-FT'' 
and replacing G^T with G'^-T - C'Of. All such operations 
can be seen to be linear transformations of the vectors of 
functions F"' and G'"'. In the appendix it is shown that 
this always leads to the same inversion procedure provided 
that the two latitudinal target functions arc equal up to a 
multiplicative constant. One implication of this is that our 
modified ]R^(g)lR^ inversion cannot be obtained from a linear 
transformation of Sekii's method, since such transformations 
never lead to a radial inversion using just F"' as Sekii's 
method docs. 

The advantage of both ours and Sekii's IR^ (g) IR^ tech- 
nique for SOLA type inversions is that it is never necessary 
to invert a matrix with a dimension of the order of the total 
number of frequency splittings (i.e. ^ 70000). The largest 
matrix to invert for construction of the latitudinal averaging 
kernels is Imax which is 250 for the mode set under consider- 
ation. For the radial averaging kernels it is in principle the 
number of different n, /-combinations which is 1380 for the 
Hare and Hounds mode set. However the number of modes 
for the radial inversions can be reduced by projection onto a 
suitable basis after performing an SVD reduction (cf. PTl; 
Christensen-Dalsgaard & Thompson 1993) to 87. 

Matrix inversion is an 0{N^) process and therefore 
naively the speed-up of the inversion due to the technique 
described here for this mode set would be a factor close to : 



modes used in the 2-D and ® R'^ inversion 
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This estimate of the speed-up was not experimentally tested. 

Finally it should be noted that a very preliminary ver- 
sion of this method was used in the GONG Hare and Hounds 
exercise (Gough & Toomre, 1993). The results there were 
poorer than those presented here, since the weighting factors 

and the scheme for obtaining the optimal resolution at 
all radii had not been satisfactorily developed at that time. 
The quality of the results is sensitive to these choices. 
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(36) 



where Nmil) is the number of observed modes with distinct 
m values for each I. However the integrations in the calcu- 
lation of the vector v in (36) turn out to dominate in the 
computing time and this is an 0{N) process so the overall 
speed up is more modest. Taking into account that the full 
2-D inversion would require computing 2-D integrals in the 
determination of v in (36) instead of 1-D integrals the the- 
oretical speed up is a factor of ~ 14 x 500 = 7000. Here 
500 is the approximate number of grid points in the radial 
direction which one would also have to integrate over in com- 
puting a complete 2-D integral (to resolve the mode kernels 
adequately) , and which can now be omitted in all 1-D latitu- 
dinal inversions. Since there is an integration for each mode 
there is an extra factor 14 from the ratio of the number of 



APPENDIX A: INVARIANCE OF KERNELS 
UNDER LINEAR TRANSFORMATIONS 

The 2D R'^ g) R^ splitting of the kernels of equation (4) is 

invariant under any linear transformation. The transformar 
tions of the vector of functions F is written as : 



fJ \i r) \F2 

and the transformation of the vector of functions G as : 



(Al) 



G2 



a b 
c d 



Gi 
G2 



(A2) 



which covers all the transformations considered in this par- 
per. The inner product of the F and G vector should not 
be affected by the transformation, therefore 



F G 



o q 
p r 



Gi 
G2 



(A3) 
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must be equal to : 

(A4) 



Fi\ Gi 



F2 J y G 

As a consequence of this it is clear that the two matrices 
must be each others inverse : 

o q\ _ 1 ( d —b 



I - ^ ^ \ I (A5) 

p r J [ad-bc) \-c a ' 

The single kernel used for the radial inversions is Fr : 

F„ = (A6) 

Here the constant (, is the ratio of the values of the kernels Q 
evaluated at the target location u = uq. In the formulation 
used in the main text 

C = -(^)' (A7) 
In the transformed set 

In this case the new radial kernel Fn becomes 
(a + C6)Ffl = ^-^l^[(a + C&)(dFi-cF2) 



+(c + Crf)(-&Fi+aF2)] 

= f 1 + C,F2 



(A9) 



The factor (a + C,b) on the left-hand side drops out again 
because of the subsequent normalization of the kernels. The 
radial part of the inversion therefore is not affected by the 
linear transformation. The latitudinal part of the inversion 
docs appear different but if the two components of the target 
function are equal up to a constant factor of multiplication 
for the original set {G\,G2), then the transformation will 
preserve that and therefore also the inversion. 
The original method as described in the paper by Sekii 
(1993b) which uses the radial kernel Fb. ~ Fi is not obtain- 
able from an invertible linear transformation of the kernels. 



